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We discuss observational consequences of f(R) dark energy scenarios that satisfy local gravity 
constraints (LGC) as well as conditions of the cosmological viability. The model we study is given 
by m(r) = C(— r — l) p (C > 0,p > 1) with m — Rf,RR/f,R and r — —Rf,ii/f, which cover viable 
f(R) models proposed so far in a high-curvature region designed to be compatible with LGC. The 
equation of state of dark energy exhibits a divergence at a redshift z c that can be as close as a few 
while satisfying sound horizon constraints of Cosmic Microwave Background (CMB). We study the 
evolution of matter density perturbations in details and place constraints on model parameters from 
the difference of spectral indices of power spectra between CMB and galaxy clustering. The models 
with p > 5 can be consistent with those observational constraints as well as LGC. We also discuss 
the evolution of perturbations in the Ricci scalar R and show that an oscillating mode (scalaron) can 
easily dominate over a matter-induced mode as we go back to the past. This violates the stability 
of cosmological solutions, thus posing a problem about how the over-production of scalarons should 
be avoided in the early universe. 



I. INTRODUCTION 

The origin of dark energy (DE) has persistently been 
one of the most serious problems in cosmology 0, Q. 
Many DE models have been proposed so far, but we have 
not found any strong evidence to support that such mod- 
els are better than cosmological constant. Thus the first 
step towards the understanding of the origin of DE is to 
find the departure from the ACDM model. 

The simplest modification to the ACDM model is per- 
haps so-called f(R) gravity in which the Lagrangian is 
written in terms of the function of a Ricci scalar R. 
It is well known that inflationary expansion is realized 
by the Starobinsky's model with a Lagrangian density 
f(R) = R + aR 2 0. Since the R 2 term is negligibly 
small relative to R at the present epoch, this model is 
not suitable to explain present accelerated expansion of 
the Universe. Instead the model with a Lagrangian den- 
sity f(R) — R— a/R n (a > 0,n > 0) was proposed to 
give rise to a late-time accelerated expansion in the met- 
ric formalism [J] (see also Refs. 0, However it was 
shown that this model is plagued by a matter instabil- 
ity @ as well as by a difficulty to satisfy local gravity 
constraints 0. Moreover it does not possess a standard 
matter-dominated epoch because of a large coupling be- 
tween dark energy and dark matter (see Refs. [H| for 
recent works). 

In Ref. [ll[ several conditions for the cosmological vi- 
ability of f(R) dark energy models were derived without 
specifying the forms of f(R). This can be well under- 
stood by considering a trajectory of each model in the 
(r,m) plane, where r = -Rf tR /f and m = Rf,Rn/f,R. 
The existence of a saddle matter-dominated epoch re- 
quires the conditions m > and — 1 < dm/dr < around 
the point (r, m) = (—1,0). The matter era can be fol- 
lowed by a stable de-Sitter attractor on the line r = — 2 
provided that < m(r = —2) < 1. This method is useful 



to rule out some of f(R) models such as f(R) = R—ot/R n 
(a > 0, n > 0) easily. 

More recently a sequence of cosmologically viable f(R) 
models was discussed in Ref. [H, EH • O ne °f sucn mod- 
els, for example, is f(R) = (R b - h) 1 ' b with < b < 1, 
which corresponds to a straight line m(r) = (b—l)(r+ 1) 
connecting the matter point (r, m) = (—1,0) to the de- 
Sitter point on the line r = — 2. The parameter m that 
characterizes the deviation from the ACDM model is con- 
strained to be m < O(0.1) for such models from the data 
of Supernova la (SN la) and Cosmic Microwave Back- 
ground (CMB) . Meanwhile local gravity experiments 
constrain the value of m to be very much smaller than 
unity in high-density regions where gravity experiments 
are carried out. This means that the deviation from the 
ACDM model needs to be very small in a high-curvature 
cosmological epoch whose Ricci scalar R is much larger 
than the present cosmological value i?o- 

A number of authors [3, EE EH recently proposed 
f{R) dark energy models that can satisfy both cosmo- 
logical and local gravity constraints (LGC) by u sing the 
so-called chameleon mechanism [l?} (see Refs. [l8l EH 
for related works). From the requirement of LGC these 
behave as close as the ACDM model during radiation 
and matter dominated epochs (i? 3> Ro)- The deviation 
from the ACDM model becomes significant after the end 
of the matter era with the growth of the quantity m. 
These models satisfy the relation f(R = 0) = 0, imply- 
ing that cosmological constant disappears in a flat space- 
time. We note, however, that the Ricci scalar is frozen 
at a value R = Ri > if the solutions are trapped by 
stable de-Sitter attractors responsible for the late-time 
acceleration. Thus in these models the system does not 
reach the region R = in an asymptotic future. 

In this paper we shall study observational consequences 
of f(R) models that satisfy LGC in addition to conditions 
of cosmological viability. The models we consider are 
given by m(r) = C(— r — l) p with C > and p > 1, 
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which cover viable models proposed in literature [lj, UM, 
[l6j in the region R S> Rq- The quantity m(r) is, for 
large p, vanishingly small during matter and radiation 
epochs (r as —1), but grows to the order of C as the 
solutions approach de-Sitter attractors on the line r = 
—2. These models exhibit peculiar evolution of the DE 
equation of state, as we will see later. Moreover matter 
density perturbations evolve differently compared to the 
ACDM cosmology for rcdshifts below a critical value Zk- 
This property can be used to place constraints on model 
parameters in addition to constraints coming from LGC, 
SN la and CMB. 

This paper is organized as follows. In Sec. HI] we 
present all conditions viable f(R) DE models need to sat- 
isfy. In addition to f(R) models studied so far, we shall 
propose another model satisfying these conditions. In 
Sec. IIIII the evolution of DE equation of state and result- 
ing observational consequences are discussed in addition 
to constraints coming from the sound horizon of CMB. 
In Sec. [TV] we study how matter perturbations evolve on 
sub-horizon scales and put constraints on model param- 
eters from the difference of spectral indices of the power 
spectra between CMB and galaxy clustering. We also 
discuss the evolution of the perturbation SR and show 
that an oscillating mode called scalaron Q easily domi- 
nates over the background value R when we go back to 
the past. This generally violates the stability condition of 
f(R) models, which gives rise to another problem about 
how to avoid the over-production of scalarons in the early 
universe. We conclude in Sec. El 



II. MODELS THAT SATISFY COSMOLOGICAL 
AND LOCAL GRAVITY CONSTRAINTS 

Let us begin with the following action 
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where k 2 = 8irG (G is a bare gravitational constant). In 
what follows we use the unit k 2 = 1, but we restore the 
gravitational constant when it is needed. Note that £ m 
and £ ra d are the Lagrangian densities of dust-like matter 
and radiation, respectively, which satisfy usual conserva- 
tion equations. In the flat Fricdmann-Robcrtson- Walker 
(FRW) background with a scale factor a the Ricci scalar 
is given by R — 6(2H 2 + H), where H = a/a is & Hubble 
parameter and a dot represents a derivative with respect 
to cosmic time t. 

There are a number of constraints viable f(R) models 
need to satisfy. First of all, to avoid anti-gravity, we 
require the condition = d//di? > 0. Modified f(R) 
models possess a scalar particle whose effective mass is 
given by 
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in the regime M 2 (R) > R [H, QJ, . In order to avoid 
that the scalaron becomes tachyons or ghosts, we require 
f t RR = d 2 //di? 2 > in this region. Note that this con- 
dition can be also derived by considering the stability of 
perturbations [20l.[2l|. 

The conditions for the cosmological viability of f(R) 
models have been studied in Ref. in great details. 
This can be well understood by considering two quanti- 
ties: 
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(3) 



The ACDM model, f(R) = R-2A, corresponds to m = 
and r = —R/ (R — 2A). The quantity m characterizes the 
deviation from the ACDM model. The cosmological via- 
bility of such models is known by plotting corresponding 
curves in the (r, m) plane. 

In what follows we shall consider cosmological evolu- 
tion that starts from a radiation epoch with large and 
positive R followed by a matter era and eventually ap- 
proaches a de-Sitter attractor with R = Ri > in fu- 
ture 1 . In cosmologically viable models we study, the 
quantity m is always smaller than 1 with fn of or- 
der unity before reaching a de-Sitter attractor. Since 
l/f,RR ^> R m such cases, one can use the scalaron mass 
given in Eq. Hence the stability conditions are given 
by EF 



f,R>0, f,RR>0, for R>Ri 



(4) 



The matter-dominated point, i\f, exists on the line 
m = —r — 1 with m close to 0, i.e., (r, m) w (—1,0). 
The presence of a viable saddle matter era demands the 
conditions [ll| (see also Ref. 



m(r « —1) > , and 



K?(rw-1)<0. (5) 
dr 



If the condition ([4]) is satisfied then the variable m is au- 
tomatically positive. The second requirement in Eq. ([5|) 
implies that m(r) curves should be be present between 
the lines m = and m = — r — 1. 

There is a stable de-Sitter fixed point that leads to a 
late-time acceleration: 



Pa 



-2, < m< 1 



(6) 



If a m(r) curve staring from Pm has an intersection 
point with a line r = —2 in the region < m < 1, 
the corresponding f(R) model is regarded as cosmolog- 
ically viable. The ACDM model is a straight line that 
links P M : (r,m) = (-1,0) with P A : (r, m) = (-2,0). 
In this paper we do not consider another accelerated 
fixed point Pg that exists on the line m = — r — 1 with 



1 During the radiation era the Ricci scalar evolves as R oc t 3 / 2 
because of the presence of non-relativistic particles. 
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(V3 — l)/2 < m < 1 [IT|. This corresponds to the case 
in which R continues to decrease in future, which can 
violate the stability condition 

A number of f(R) models satisf ying the above condi- 
tions were considered in Refs. [H, [H|. Some examples 
are 



(i) f(R) = (R b 

(ii) f(R) = R- 



-A) c 
aR n 



(c>l, be ml), (7) 
(a>0, 0<n<l), (8) 



c)/c]r 



1 and 



which correspond to m(r) = [(1 
m(r) = n(l + r)/r, respectively. 

Let us next consider local gravity constraints on f(R) 
dark energy models. The LGC are satisfied for Ml 3> 1 
[T3I . f]~8l | , where I is a scale at which gravity experiments 
are carried out. Using Eqs. ^ and ©, this constraint 
is expressed by 



(9) 



where i? s is a curvature measured on the local structure 
and is proportional to the energy density p s of the struc- 
ture (R s f» 87rGyO s ). Using the present cosmological den- 
sity po and the Hubble radius H^ 1 ~ 10 28 cm (in what 
follows we use the subscript "0" for present values), the 
above constraint is rewritten as 



m(R s ) « BL 



(10) 



where we used f r b ~ 1 and R ~ fl^ ^ %<xGpQ. The 
r.h.s. of Eq. (jTUJ) is very much smaller than unity fl3T | 
because I <C -ff^ 1 even though p s is larger than p . In 
the case of the Cavendish-type experiments the typical 
constraint is m(R s ) <C 10 -43 , as we will see later. 

The above argument shows that in the high-curvature 
region (R S> Ro) the quantity m needs to be negligi- 
bly small. Cosmologically this means that during ra- 
diation and matter eras the models need to mimic the 
ACDM model with a high-precision. Note that the mod- 
els (i) and (ii) given in Eqs. ||7J) and © behave as 
m(r) = C(—r — 1) as r approaches — 1. In such cases, 
however, LGC are not satisfied unless C is chosen to be 
unnaturally small. 

Hu and Sawicki [l4| proposed an explicit f(R) model 
that satisfies both cosmological and local gravity con- 
straints. It is given by 



f(R) =R- XR C 



(R/Rc) 
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where the power 2n is used instead of n. Starobinsky [If 
also proposed another viable model: 



f(R) =R- XR C 



1 R2 

1 + m 



(12) 



In both models n, A and R c arc positive constants, where 
R c is the order of the present Ricci scalar Rq. Since 



f(R = 0) = cosmological constant disappears in a flat 
spacetime. Thus the origin of dark energy can be re- 
garded as the geometrical one. 

Let us check the cosmological viability as well as the 
stability for such models. In the region R R c these 
behave as 



f(R) 



R — XR r 



R 

2n(2n+ 1) 
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(13) 

(14) 
(15) 



Thus in this region the models (fTTj) and (fT2"]) have the 
following property 



m(r) = C(-r - If 



(16) 



where p = 2n + 1 > 1 and C is a positive constant. It 
is obvious that, for larger p, m(r) becomes very small as 
r — > — 1 so that the model satisfies LGC. Since ^(r = 
— 1) = the condition ([5]) is also satisfied. 

Let us next check the conditions ^ and ©. In the 
model (fTTj) , the de-Sitter point at r = —2 is determined 
by the value of A: 



[l + x\ n ) 2 



„1n- 



i(2- 



2a; 2 ™ 
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(17) 



where x\ = R\/R c . From the stability condition < 
m(r = —2) < 1 wc obtain 

2x{ n - (2n - l)(2n + <±)x\ n + (2n - l)(2n - 2) > .(18) 

When n = 1, for example, we have x\ > \[2> and 
A > 8V3/9. Under Eq. (Tjj> one can show that the con- 
dition is satisfied. This situation is similar to the 
Starobinsky 's model (fT2"j) . see Ref. [lfl] for details. 

We can extend the above two models to the more gen- 
eral form 

f(R)=R-£(R), £(0) = 0, £(i? > Rc) -> const. (19) 

The conditions (j4]) translate into 

Cr<1, iRR<0, for R>Rl (20) 

In order to satisfy LGC, we require that £(i?) approaches 
a constant rapidly as R grows in the region R 3> R c 
(such as £,(R) — constant — (R c /R) 2n discussed above). 
Another model to meet these requirements is 



f(R) = R- Ai? c tanh ( — 
Rc 



(21) 



where A and R c are positive constants. A similar model 
was proposed by Appleby and Battye [l6| , although it is 
different from (|2T|) in the sense that £,(R) can be negative 
for R < Ri. In the region R > R c the model {T3TJ) 
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behaves as f(R) ~ R— Ai? c [l — exp(— 2R/R C )], which can 
be regarded as the special case of (fTrJ)) with a limit p — > 
oo. The Ricci scalar at the de-Sitter point is determined 
by A, as 

A= , ^ cosh2 (^) , (22) 

2 sinh(a;i) cosh(a;i) — x\ 

where x\ = R\/R c . From the stability condition ^ of 
the de-Sitter point we obtain the constraint 

xi> 0.920, A > 0.905. (23) 

In the models (fTT|) and (fT2|) . f^RR are negative for 

< R/R c < [(2n - l)/(2n + l)] 1 / 2 " and < R/R c < 

1 / y / 2n + 1, respectively. In the model pip the quantity 
f,R also becomes negative (i.e., > 1) in some regions, 
whereas in the model (fT2"|) it is possible to have f r > 
for all positive R in some restricted regions of the pa- 
rameter space (0.944 < A < 0.966 for n = 2). See the 
curves (i) and (ii) in Fig. Q] for illustration. Of course 
we are considering the situation in which the violation of 
the conditions f t R > and f y RR > occurs for R < R\, 
so it is harmless as long as the universe is in the region 
R> R\. In the model (|2"Tj) we always have J,rr > for 
positive R, whereas f t R is positive for A < 1. Hence, if 
0.905 < A < 1, this model satisfies the conditions f t R > 
and f t RR > for all positive R and also possesses a de- 
Sitter attractor. Such an example is plotted as the case 
(iii) in Fig. [2 

In Fig. [5] we plot the trajectories in the (r, m) plane 
for the model ([12]) with n = 2, A = 0.95, and for the 
model (|21[) with A = 0.95. The solutions start from the 
region R/R c 3> 1 around the point Pm and they finally 
approach the de-Sitter point Pa &t R = R\. Since f t R > 
(R > 0) and /(0) = in such cases we have f(R) > 
for positive R, which means that r = —Rf,R/f is always 
negative. The quantity m = Rf .rr/ f ,R is positive for 
R > Ri because /,rr > 0. In the case (A) of Fig. [3] 
we have m < for < R/R c < \j\J2n + 1 because f t RR 
changes the sign. Meanwhile, in the case (B), m is always 
positive for R > 0. In the limit R/R c — > the trajectories 
approach the point (r, m) = (—1,0) again in both models. 
If the Ricci scalar oscillates around R = 0, the quantity 
f,RR becomes negative for R < even for the model (f2"Tj) . 
As we will see later in details, this can indeed occur by 
the oscillation of scalarons unless initial conditions are 
appropriately chosen. The above argument shows the 
importance of confining the Ricci scalar in the region 
R > Ri to ensure the stability of models. 

In this paper wc shall study a number of cosmological 
constraints on the models of the type (flB"]) . As we men- 
tioned, in the region R ^> R c , this covers the models (fTTj) 
and (TT2"]) with p = 2n + 1 as well as the model (|2ip with 
the limit p — > oo. If < C < 1 there exists a stable de- 
Sitter point P A at r = -2. The model m(r) = C(-r-l) p 
essentially contains sufficient information about how vi- 
able f(R) models behave. During radiation and matter 
eras (r ~ —1) the quantity m is very much smaller than 
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Figure 1: The illustration of £(i?) as a function of R/R c in 
three different models. Each corresponds to (i) the model 
by Hu and Sawicki with n = 1 and A = 2, (ii) the model 
(|12p by Starobinsky with n = 2 and A = 0.95, and (iii) the 
model ([211 with A = 0.98. We also plot a line £(i?) = R to see 
whether or not the condition < 1 (i.e., > 0) is violated. 
The black points represent de-Sitter fixed points (R = Ri). 
Note that in the case (i) the de-Sitter point corresponds to 
Ri/Rc = 3.383, which is outside of the figure. In the case (i) 
we have f.R < for 0.296 < R/R c < 1, whereas in the cases 
(ii) and (iii) f t n > for all positive R. In the models (JTTJ) 
and (|12p it is inevitable to avoid that becomes negative 
in the small R region, but in the model pip it is possible to 
realize J,rr > for R > 0. In the region before the solutions 
reach the de-Sitter attractor (i.e., R > Ri), both f t n and /,rr 
are positive in the above three models. 

unity, but it grows to the order of C once the system 
approaches the de-Sitter attractor Pa- Thus one can see 
the departure from the ACDM model around the present 
epoch. Note that we are only concerned with the region 
of m(r) curves before the solutions reach the de-Sitter at- 
tractor. If we demand the condition f(R = 0) = 0, this 
can be satisfied by modifying the form of m(r) outside 
the region that connects Pm to Pa (as in the models (TT2"]) 
and (HU in Fig. HJ). 

Before entering the details of various cosmological con- 
straints, we consider LGC on the model (TTu| . In high- 
dense regions where local gravity experiments are car- 
ried out (R s ^S> R c ~ Rq) the quantity r behaves as 
Eq. (fnj and hence m ~ rh(Ro/R s ) p , where m is a con- 
stant whose order is not much different from unity. Since 
Rq/R s ~ Po/Ps, the constraint (fit)]) yields 

In the Cavendish-type experiments the typical values are 
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Figure 2: Two trajectories in the (r, m) plane. The trajectory 
(A) corresponds to the model (|12[) by Starobinsky with n — 2 
and A = 0.95, whereas the trajectory (B) to the model (|21[) 
A = 0.95. In both cases the solutions start from the region 
around the point Pm '■ (r,m) — (—1,0) with R ^> R c . They 
approach the stable de-Sitter point Pa on the line r = —2 
with < m < 1. In the case (A) the quantity f : RR becomes 
negative for < R/R c < l/y/2n + 1, while f.R is positive for 
R > 0. In the case (B) one has /,rr > and f t R > for all 
positive R. In the limit R/R c — > both models approach the 
point (r,m) = (—1,0) again. 



g/cm 3 and 



10- 



cm ]2c 

-l 



Recalling the 



Ps ~ 10" 

values po ~ 10~ 29 g/cm 3 and iT 1 ~ 10 28 cm, we find 
that the LGC is well satisfied for p > 3. 

In the case of solar-system experiments, if we take the 
typical solar-system length scale I = 1 Au~ 1.5 x 10 13 cm 
with the density p s ~ 10 -24 g/cm 3 at the distance from 
Sun [l4||, we obtain the constraint p > 5. Meanwhile 



the so-called Shapiro time-delay effect |24| conies mainly 
from the gravity contribution around the radius of Sun 
(£ ~ 7.0 x 10 10 cm) with the density p s ~ 10~ 15 g/cm 3 , 
which gives a much weaker constraint: p > 2. Thus the 
constraint p > 5 is certainly enough to satisfy LGC and 
is even too tight in some of gravity experiments. 



III. SN IA AND CMB SOUND HORIZON 
CONSTRAINTS 



In this section we discuss the cosmological evolution 
of the model ([To]) at the background level and confront 
it with constraints coming from SN la and the sound 
horizon of CMB. In the flat FRW spacetime the variation 
of the action (JTJ) leads to the following equations 



ZFH" 



Prad + (FR - f)/2 - 3HF, (25) 



-2FH = p m +(4,/3)p rad +F-HF, 



(26) 



where F = df /OR. Here p m and p ra d are the energy den- 
sities of a non-relativistic matter and radiation, respec- 
tively, which satisfy the usual conservation equations. 

Following Refs. [ijlGJ] we introduce the dimcnsionless 
variables 



F 
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Xl HF ' X2 6FH 2 



X3 



R 



Q H 2 ' X4 3FH 2 ■ 

Then we obtain the dynamical equations (Til 
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— 1 — X3 — 3^2 + X 
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x 2 {2x 3 - 4 - xi) , 



2a; 3 (a;3 ~ 2) 



m 

-2x3X4 + X1X4. 



(27) 

(28) 
(29) 

(30) 
(31) 



where a prime represents a derivative with respect to 
TV = In a. Since m is a function of r = x 3 /x2, the above 
system is closed. The energy density and the pressure 
of DE to confront with SN la observations are given in 
Ref. [l3[ and the corresponding equation of state (EOS) 
of DE is 



1 



2x 3 - 1 + (F/F ) Xi 



31-{F/F ){l- Xl -x 2 



Xi) ' 



(32) 



where Fq is the present value. 

Our model (p~6|) needs to satisfy the condition > 
for R > Ri , which leads to the increase of F toward the 
past as R gets larger. The denominator in Eq. (|32|) is 
written as 1 — (F/ Fo)il m , where 



Pn 



3FH 2 



1 — Xl 



X2 — X3 — Xi . 



(33) 



Since f2 m increases from present to the matter-dominated 
epoch, it happens that u>de exhibits a divergence at a 
redshift z c satisfying f2 m = Fq/F. This is in fact generic 
to cosmologically viable models that fulfill the criterion 
(|4]). The EOS of DE crosses a cosmological constant 
boundary (wde = —1) at a redshift Zb smaller than z c 
0. 

In Fig. [3] we plot u>de versus the redshift z = a$/a — 1 
for the model m(r) = C{— r — l) 2 with three different 
values of C . In all simulations the present epoch (z = 0) 
is identified as a matter energy fraction Q m = 0.28 with a 
radiation contribution f2 ra( j = p ra d/3FiJ 2 ~ 10~ 4 . Note 
that the current universe is in the middle of approach- 
ing the de-Sitter attractor Pa- (xi,X2,X3) = (0,-1,2) 
from the matter point Pm- (xi,X2, X3) ~ (0, —1/2, 1/2). 
As we see in Fig. [31 wde is larger than — 1 at z = 
and decreases to —00 as z — > z c — after crossing the 
cosmological constant boundary at z = z/,. For smaller 
C, z c gets larger. When C — 1 we have z c = 2.03 
and wd~e(z = 0) = —0.853, whereas if C = 0.5 we get 
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Figure 3: The DE equation of state wde versus the redshift z 
for the model m(r) = C( — r — l) 2 with three different parame- 
ters (C = O.f , 0.5, f .0). With the increase of C the divergence 
of ibde occurs for smaller z. 



z c = 2.31 with wde(z = 0) = -0.904. This peculiar 
behavior of the EOS of DE is an interesting signature 
to discriminate f(R) models from the ACDM cosmology. 
In particular the larger deviation from the ACDM model 
leads to a smaller critical redshift z c that can be reached 
in future observations. 

For the values of p greater than 2, z c gets larger and 
wde(z = 0) tends to be closer to —1. In Table I we show 
Zb, z c and u>de(z = 0) together with present values of 
to for several different choices of p and C. When p > 5, 
if we look at the low-redshift region only, the models 
are hardly distinguishable from the ACDM model. Still 
the EOS of DE shows a peculiar behavior in the high- 
redshift region: z > 3. For larger p the cosmological 
constant boundary crossing occurs at the redshift close 
to the present epoch. We note that the recent SN la data 
analysis finds some evidence for such a crossing p5j . If 
future high-precision observations favour models whose 
EOS corresponds to a phantom (wde < — 1) in most of 
the past epochs relevant to SN la observations, this can 
be the signal of f(R) gravity. 

If we use the criterion w-£>e(z = 0) < —0.7 accord- 
ing to the current SN la data (26|, we find from Table 
I that the models with p > 2 satisfy this requirement 
even when C is as close as unity. The slopes of the EOS, 
|du>DE/dz|, are found to be smaller than the order of 0.1 
around the present epoch, which do not provide addi- 
tional information to constrain models. When C = 1, 
we have wde(.z = 0) > —0.7 only for p < 1.4. Thus the 
current SN la observations do not provide a better con- 
straint on the power p than the one obtained by LGC. 



p 


C 




Zc 


ihde(z = 0) 


m{z = 0) 


1.5 


0.1 


0.79 


3.44 


-0.959 


0.075 


1.5 


0.5 


0.73 


2.24 


-0.851 


0.244 


1.5 


1 


0.85 


2.14 


-0.743 


0.293 


2 


0.1 


0.62 


3.88 


-0.969 


0.067 


2 


0.5 


0.54 


2.31 


-0.904 


0.220 


2 


1 


0.56 


2.03 


-0.853 


0.285 


3 


0.1 


0.40 


5.00 


-0.982 


0.055 


3 


0.5 


0.33 


2.83 


—0.955 


0.180 


3 


1 


0.32 


2.31 


-0.936 


0.248 


5 


0.1 


0.16 


7.25 


-0.994 


0.037 


5 


0.5 


0.10 


4.30 


-0.990 


0.124 


5 


1 


0.09 


3.53 


-0.988 


0.180 



Table I: The values Zb, z t 
model m(r) = C(— r — 1 
and C. 



, woe{z = 0) and m(z = 0) in the 
| p for several different choices of p 



It will be interesting, however, to carry out a likelihood 
analysis using the future data of SN la along the line of 
Refs. 

Let us also consider the sound horizon constraint com- 
ing from CMB. The angular size of the sound horizon is 
defined by 



3 s (z)dz 
H(z) 




dz 
Wz) 



(34) 



where c 2 s (z) = l/[3(l + 3/0&/4p^)] is the adiabatic baryon- 
photon sound speed and Zdcc — 1089. From the po- 
sition of CMB acoustic peaks we obtain the constraint 
6 S = 0.5946 ± 0.0021 deg from the WMAP 3-ycar data 
[27| . For the models in which the effect of dark energy 
is not negligible during the matter-dominated epoch, the 
quantity O s is rather strongly modified (as in the coupled 
quintessence [28|). 

In f(R) gravity, if the quantity to is not much smaller 
than 1 during the matter era, this leads to a considerable 
change of O s compared to the ACDM model. In fact 
this happens for the models ([7]) and ([5]), as was shown 
in Ref. [13J. In our model (| 16[) the quantity to is very 
much smaller than unity during the matter era from the 
requirement to satisfy LGC. Hence it is easier to satisfy 
the sound horizon constraint compared to the models ([7]) 
and ([8]). In fact we have evaluated Q s numerically and 
confirmed that the models with p > 2 are consistent with 
the WMAP 3-year data. Thus the data coming from the 
CMB sound horizon does not provide tighter constraints 
relative to the SN la data. 



IV. MATTER PERTURBATIONS AND 
SCALARON OSCILLATIONS 



In this section we study constraints on the model |TJ 
coming from matter density perturbations. Let us con- 
sider scalar metric perturbations a, /3, tp and 7 about the 
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fiat FRW background [H: 

ds 2 = -(1 + 2a)ctt 2 - 2a(3 ii dtdx i 

+a 2 [(1 + 2<p)<5y + 2 7 |y] darW . (35) 

In what follows we neglect the contribution of radiation 
as it is unimportant to discuss the evolution of matter 
perturbations during the matter-dominated epoch. The 
energy-momentum tensors of a pressureless matter is de- 



composed by 

Tq = -(p m + SPm) , T? = -p m V m)i , (36) 

where v m is a velocity potential. 

Introducing a covariant velocity perturbation, v = 
av m , we obtain the following equations of motion in the 
Fourier space [3l[ (see also Refs. pH I32I]): 



a = v , 
(5p m /p m ) 

k + 2Hk ■ 



k 2 

k — 3Ha -v . 

















2F 





k 2 



R 



1 



1 



SF + 3HSF +1— — — \5F= -Sp m + F(n + a) + (2F + 3HF)a - -FSR 



X + H X -a-<p= ^(5F~F X ), 



3H5F 



k 2 

3H + 3H 2 ) 5F - 3HFa -Fk- 5p r , 



(37) 
(38) 

,(39) 

(40) 

(41) 

(42) 



r 



where k is a comoving wavenumber and K = i{Ha — 
ip) + ([3 + (i'y)k 2 /a. We define a gauge-invariant quantity: 
8 m = 6p m /p m + 3Hv. In the comoving gauge where 



v = 0, wc find from Eqs. (|37[) and (f3"8")l that a = and 
k = S m . Then from Eqs. (|3"9"|) and (|^U|) we obtain 



SF + 3HSF + ( — 



Pm x 1 
:0„, = 



2.F 1 



2F 



-6H 2 



3F, 



- 4iP - 2H )6F = -5p 



1 



5F + 3H5F + 3SF 



F5„ 



(43) 
(44) 



In the model (fT6)) the quantity m = R/^rr/F is very 
much smaller than unity during the matter era with 
F ~ 1. Since 1/ / jm 3> i? in such the scalaron 

mass squared is given by Eq. © and satisfies the rela- 
tion M 2 ^> R ~ i7 2 . In what follows we shall discuss 
two cases: (A) M 2 > k 2 /a 2 and (B) M 2 < fc 2 /a 2 , sepa- 
rately. As we will see below, the modes that are initially 
in the region (A) can enter the region (B) during the 
matter-dominated epoch. 



A. The region M 2 > k 2 /a 2 

When M 2 k 2 /a 2 , Eq. (f44|) is approximately given 
by 

SF + 3HSF + M 2 SF ~ \s Pm , (45) 

where we used the fact that the variation of the quantity 
F is negligibly small during the matter era. This is a 
very good approximation for the model (|16[) . since m is 
vanishingly small during the matter era. 

The general solutions for Eq. (|45[) are given by the 
sum of the oscillating solution SF osc obtained by set- 
ting Sp m = and the special solution SFi n d of Eq. (|4"5"|) 
induced by the presence of matter perturbations Sp m . 
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The former was obtained by Starobinsky [l5j for the 
model (fT2")) in the unperturbed flat FRW background 
(i.e., k = 0). The oscillating part SF osc satisfies the 
equation (a 3/2 <5F osc )" + M 2 (a 3 / 2 5F osc ) ~ 0. By using 
the WKB approximation, we obtain the solution 

<5F osc oc a' 3 / 2 fW' 4 cos ( J ^|== d ^ ■ (46) 

During the matter era in which the background Ricci 
scalar evolves as = 4/(3i 2 ), the quantity f^M. has 
a dependence f.RR oc oc t 2 ^ v+1 \ Hence the evo- 

lution of the perturbation, SR osc = 5F osc / f rr, is given 
by 

SR osc ~ct — cos(c i p ), (47) 

where c and cq are constants. As we go back to the past 
the amplitude of <5i? sc dominates over R^°\ unless the 
coefficient c is chosen to be very small. Since R gets 
smaller than R\ and even becomes negative, the stability 
condition ((4| is violated. This property also holds in 
the radiation era during which 5R osc and the background 
Ricci scalar i?' ' evolve as 

SR osc ~ ct- 3 ^ costcoi-^- 1 ') , i? (0) oc r 3 ' 2 .(48) 

Thus we need to avoid the excessive production of 
scalarons in the early universe so that \SR osc \ < 
is satisfied at all times. This problem is even severe for 
the models of the type (j2"T]) . Moreover the scalaron mass 
rapidly grows to the past in these models and can exceed 
the Planck mass even during the matter era. 

The special solution SF ind of Eq. (|4"5)l can be derived 
by using the approximation used in Refs. [H, [Hj]. This 
amounts to neglecting the first and second terms relative 
to others, giving 

SFi n d ~ f.RR fiPm , SRind ~ 5p m ■ (49) 

Under the condition |£F osc | <C |<5-Fi n d| we have 5F ~ 
f.RR&Pm- Substituting this relation for Eq. (|43|) and us- 
ing the property M 2 3> k 2 /a 2 , we obtain 

<L + 2H6 m - AnGprnSm ^ . (50) 

Here we have reproduced the gravitational constant for 
clarity. This is the usual equation of matter perturba- 
tions on sub-horizon modes in ACDM cosmology and has 
a growing mode solution S m oc a oc t 2 / 3 . From Eq. (J49J) 
we get 

SF ind cx t 2 P+ 2 ' 3 , 8R ind ^t~ 4 / 3 . (51) 

Compared to the oscillating mode (|47| , the induced mat- 
ter mode 5Ri n d decreases more slowly and thus dominates 
in the late universe. Relative to the background value 
R(°>, the perturbation, 5R = 5R OBC + <5i?i n d, evolves as 

~ h t~ ^ cosfoi-") + b 2 1 2 ' 3 , (52) 



T 



o.oio r ~ m 




5 10 15 20 25 30 35 

Z 

Figure 4: The evolution of 5R, 5R/R {0) and S m for the 
model m(r) — (— r — l) 3 with the mode k/aoHo = 335 in 
the case where the coefficient bi in Eq. (|52p is very small 
so that the scalaron mode SR osc is negligible relative to the 
matter induced mode SR[ n d. The transition from the region 
M 2 ^> k 2 /a 2 to the region Af 2 <C k 2 /a 2 occurs around the 
redshift Zk — 5. 



where 61 and 62 are constants. Unless the coefficient b\ is 
very small, the oscillating mode dominates over the mat- 
ter induced mode to violate the condition R > R± as we 
go back to the past. In Fig. 2] we show an example about 
the evolution of perturbations in which the initial condi- 
tion of SR is chosen to be very close to Sp m (see Appendix 
for perturbation equations suitable for numerical calcu- 
lations). The perturbation evolves as SR oc i~ 4 / 3 during 
the period in which the condition M 2 3> k 2 /a 2 is satis- 
fied. After the system enters the region M 2 <C k 2 /a 2 , SR 
decreases more rapidly as we see in the next subsection. 
Since SR ~ Sp m and R^ ~ SH 2 ~ p m during the matter 
era, we obtain the relation SR/ R^ ~ <5 m . This property 
is in fact confirmed in Fig. 0]in the region M 2 3> k 2 /a 2 . 

Figure [5] is the case in which the oscillating mode dom- 
inates over SRi n d around the redshift z > 30. Since \SR\ 
grows to the order of R^ the Ricci scalar R becomes 
negative in this region, thus violating the stability con- 
dition These results confirm that the coefficient 61 
should be chosen to be very small to avoid the dominance 
of the scalaron mode in an early epoch. 

B. The region M 2 < k 2 /a 2 

Since the scalaron mass decreases as M oc t~( p+1 ' , the 
modes which initially exist in the region M 2 3> k 2 /a 2 
can enter the regime M 2 <C k 2 /a 2 during the matter- 
dominated epoch. In this regime Eq. (|44[) is approxi- 
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Figure 5: The evolution of 5R, SR/R^ ' and S m for the model 
m(r) = (— r — l) 3 with the mode fc/ao-ffo = 315 in the case 
where the coefficient bi in Eq. (|52p is not chosen to be very 
small. The scalaron mode 5R osc dominates over the matter 
induced mode 8Ri n d around the redshift z > 30. 



matcly given by 



V 


ZA 


Zk 


8n w (t A ) 




Sri*- N) (t ) 


2 


0.95 


9.62 


0.106 


0.107 


0.108 


3 


0.86 


4.83 


0.075 


0.074 


0.077 


4 


0.81 


3.25 


0.057 


0.056 


0.061 


5 


0.78 


2.49 


0.047 


0.045 


0.053 


6 


0.76 


2.03 


0.039 


0.035 


0.044 


7 


0.75 


1.72 


0.034 


0.028 


0.039 



Table II: The redshifts za and Zk in the model m(r) = (— r — 
l) p for the mode k/aoHo = 300. We also show analytic and 
numerical values of 5n(tA), which are denoted as Su^^a) 
and 8n^ N \tA) respectively. The values Sn^ N \to) are obtained 
by numerically integrating perturbation equations up to the 
present epoch. 



Then the evolution of the perturbation SR 
6 Ri n d, relative to R(°\ is given by 



SRo 



SR 



i hi- 2 *- 2 ' 3 cos(c fct 1 / 3 ) + b 2 t- 2 P+^ 



(59) 



As long as the scalaron mode is suppressed at the be- 
ginning of the matter era, the second term on the r.h.s. 
of Eq. (|59|) dominates over the first one. In Figs. 0] and 
[5] the sudden decrease of SR means that the system en- 
ters the region M 2 <C k 2 /a 2 in which the evolution of 

SR is characterized by SR oc i _2p+ ^~B . At this stage 
SR/R 1 - ) is no longer proportional to S m . 



SF + 3HSF +^rSF~-Sp n 



(53) 



C. The matter power spectra 



Using the WKB approximation, the solution correspond- 
ing to the scalaron mode is 



SR n 



SF 



f.RB 



^ ~ci- 2p - 8 / 3 cos(c fci 1/3 ). (54) 



The matter-induced special solution of Eq. (|53")) is ap- 
proximately given by 



(55) 



From Eq. (|43[) we obtain the following approximate 
equation under the condition |(5.Fosc| <C |&Find| : 



Sm + 2H5 m - - ■ AnGprnSm ~ . 



(56) 



Relative to the region M 2 k 2 /a 2 the growth rate of 
S m is enhanced and is given by 



^33-1 
0m. OC t 6 



Hence the induced-matter mode evolves as 



(57) 



%/33-S n.. i V33-1T 

5F ind oc. ts— , SR ind (x t- 2p+ ^5— . (58) 



The evolution of the matter perturbation is given by 
S m oc t 2 / 3 for M 2 > k 2 /a 2 and S m cx fC^- 1 )/^ for 
M 2 <C k 2 /a 2 . We shall use the subscript "fc" for the 
quantities at which k is equal to aM, whereas the sub- 
script "A" is used at which the accelerated expansion 
starts (a = 0). While the redshift za is independent of 
k, Zk depend on k and also on the mass M. 

In Table II we show numerical values of za and Zk for 
the mode k/aoHg = 300 in the model m(r) = (—r — l) p . 
For smaller p the period of a non-standard evolution of 
8 m becomes longer because Zk tends to be larger. We also 
note that, for larger k, Zk gets larger. This means that 
the duration of the period of an additional amplification 
of S m is different depending on the mode k, see Fig. [5] 
Since the time tk has a dependence oc k~ 3 p+! , the 
matter power spectrum P$ m = (k 3 /2ir 2 )\5 m \ 2 at the time 
t\ shows a difference compared to the case of the ACDM 
model: 



A 



ACDM 



(*A) 



' \/33-l 



n/33-5 
OC k 3 P + 1 



(60) 



While the galaxy matter power spectrum is modified 
by this effect, the CMB spectrum is hardly affected 
except for low multipoles around which an integrated 
Sachs- Wolfe (ISW) effect becomes important 35(. Thus 



10 




i o 

z 



20 



0.75 



10 



15 



Figure 6: The evolution of the matter perturbation <5 m in the 
model m(r) = (— r — l) 2 for the modes k/a^Ho = 6, 300, 3000. 
The transition redshift z k increases for larger k. For the mode 
k/aoHo = 300 we have z% = 9.62. 



there is a difference for the spectral indices of two power 
spectra, i.e., 



3p+l 



(61) 



Since z k becomes as close as z\ for larger p, it is not 
necessarily guaranteed that Eq. (foTj) is valid in such cases. 
Moreover the estimation (j6"T]) does not take into account 
the evolution of S m after z = z\ to the present epoch {z = 
0) . In order to see the validity of the formula (|61| we have 
evaluated numerical values of 6n(t\) as well as Sn(to) 
integrated up to the present epoch. From Table II we find 
that the estimation (|61|) agrees well with the numerical 
obtained Sn(t\) for p < 5. The difference appears for 
p > 6, but it is not significant. 

After the system enters the epoch of an accelerated ex- 
pansion, the momentum k can again become smaller than 
aM. Hence the fc-dependence is not necessarily negligi- 
ble even for z < z\. However we find from Table II 
that 5n(to) is not much different from Sn(t\) derived by 
Eq. (f6Tj) . Thus the analytic estimation (j6Tj) is certainly 
reliable to place constraints on model parameters except 
for p> 1. 

Obscrvationally we do not find any strong signature for 
the difference of slopes of the spectra of LSS and CMB 
|36j ] . If we take the mild bound Sn < 0.05, we obtain the 
constraint p > 5. It is interesting to recall that LGC are 
well satisfied for p > 5. 

Finally we shall discuss the integrated Sachs-Wolfe 
(ISW) effect in the model m(r) = C{—r — l) p . In or- 
der to confront the model with CMB it is convenient to 



Figure 7: The evolution of the effective gravitational potential 
$ = (a — <p)/2 in the model m(r) = (— r — l) p for the mode 
fc/ao-ffo = 4 with three different values of p. Note that $ 
is initially normalized as unity. We also show the evolution 
of $ in the ACDM model. For larger p the decay of the 
gravitational potential is more significant. When p > 2 the 
models are well consistent with CMB low-multipole data. 



study the evolution of an effective gravitational poten- 
tial $ = (a — tp)/2 in the longitudinal gauge (x = 0) 
[2H . |32|. Under the sub-horizon approximation used in 
RefsH33l H| we obtain, from Eqs. (gT]) and (gSJ), the 
following relation 



(62) 



In the ACDM model the gravitational potential remains 
constant during the standard matter era, but it decays 
after the system enters the accelerated epoch. This leads 
to the ISW effect for low multipoles of the CMB power 
spectrum. In our f(R) model (|16[) . the additional growth 
of matter perturbations in the region z < z k changes the 
evolution of <&. 

In Fig. [7] we plot the evolution of $ in the models 
m(r) = (— r — l) p with the mode k/aoHo = 4 for several 
different values of p. For smaller p the gravitational po- 
tential does not decay much because the duration of the 
region z < z k gets longer. The models that cancel the 
ISW effect by the additional growth of 5 m are consistent 
with the CMB data [35|, which means that the models 
with p > 2 are allowed. Hence the information coming 
from the ISW effect of CMB does not provide strong con- 
straints on the model parameters compared to the galaxy 
spectrum discussed above. 
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V. CONCLUSIONS 

We have discussed various observational signatures of 
f(R) dark energy scenarios that satisfy both cosmological 
and local gravity constraints. The f(R) models do not 
have much freedom to fulfill all such constraints. Gen- 
erally the models need to mimic the ACDM model in 
a high-curvature region where local gravity experiments 
are carried out (R 3> Ro ~ Hq)- The deviation from 
the ACDM cosmology becomes important after the end 
of the matter-dominated epoch in such viable models. 

The models given in Eqs. (|11[) and (TT2]) belong to such 
viable classes that possess a de-Sitter attractor at R = 
Ri > and satisfy the condition f(R = 0) = 0. The 
stability condition requires that f r > and f.RR > 
for R> Rl In the models (JTT]) and JT2]) either of these 
conditions is violated in the region R < i?i , but this is not 
problematic as long as the Ricci scalar does not oscillate. 
In Sec. [TTl we discussed general properties about viable 
f(R) models and proposed another simple model (J5TJ) 
that can satisfy the conditions f.R > and J,rr > for 
all positive R. 

The model we have studied in this paper is given by 
m(r) = C(— r — l) p (C > 0,p > 1), where m and r are de- 
fined by Eq. §3§ . In the high-curvature region character- 
ized by R Rq, these recover the models (TTTI) . (|T2"j) with 
finite p and also the model (j2"Tj) with p — > oo. Moreover 
the departure from the ACDM cosmology can be cap- 
tured by the growth of the quantity m as the solutions 
approach the de-Sitter fixed point on the line r = — 2. 
In fact this model possesses rich observational signatures 
relevant to SN la, galaxy clustering and CMB. 

The equation of state of dark energy shows a peculiar 
divergent behavior with a cosmological boundary cross- 
ing. When the deviation from the ACDM model is signif- 
icant, the redshift z c at which such a divergence occurs 
can be as close as a few. This may be detectable in future 
observations of SN la and weak lensing. Using observa- 
tional bounds on the equation of state of dark energy at 
present [z = 0), we found that the models with p > 2 are 
consistent with the data. The deviation parameter m is 
constrained to be m(z = 0) < 0.3. We also showed that 
the models with p > 2 satisfy sound horizon constraints 
of CMB. 

We have discussed the evolution of matter density per- 
turbations 8 m together with the perturbation in the Ricci 
scalar R. The mass squared of the scalaron, M 2 ~ 
1 / (3/,_r_r), is much larger than H 2 and can cross the value 
k 2 /a 2 during the matter era. In the early epoch with 
M 2 ^> k 2 /a 2 the matter perturbation evolves as in the 
standard way, provided that the oscillating mode 5R osc 
(scalaron) is suppressed relative to the induced matter 
mode <5i?md- However the scalaron dominates at early 
epochs unless the coefficient of this mode is chosen to be 
very small so that SR is always as close as <5i?i n d- The 
dominance of the scalaron means the violation of the con- 
dition which leads to an instability of solutions in the 
matter-dominated epoch. This property persists in the 



radiation era and hence it poses a serious problem about 
how an over-production of the scalaron is avoided in the 
early universe. 

In the late epoch characterized by M 2 <C k 2 /a 2 the 
matter perturbation evolves in a non-standard way, see 
Eq. ([56]) . This leads to additional growth of matter per- 
turbations depending on the wavenumber k. Considering 
the evolution of 5 m by the time t\ at which an accelerated 
expansion sets in, the difference about spectral indices of 
the power spectra between galaxy clustering and CMB 
is given by Eq. (frJTj). We also integrated perturbation 
equations numerically by the present epoch and found 
that the estimation (|61[) agrees fairly well with numeri- 
cal values. Using the rather mild criterion 6n < 0.05, the 
constraint on the parameter p is given by p > 5. We have 
also studied the evolution of an effective gravitational po- 
tential and found that the integrated Sachs- Wolfe effect 
in low multipoles of CMB does not provide a stronger 
constraint than the one coming from LSS. 

It will be certainly of interest to place stringent con- 
straints on the model parameters using future high- 
precision observational data. We hope that this allows 
us to find some signatures about the deviation from the 
ACDM model. In particular the detection of unusual 
behavior of the equation of state of dark energy can be 
strong evidence for f(R) gravity models. 
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Appendix 

In this Appendix we present the perturbation equa- 
tions convenient for numerical purpose. We neglect the 
contribution of radiation, i.e., x± = 0. Using the dimen- 
sionlcss variables given in Eq. Q27|) and introducing a new 
quantity SF = 6F/F, Eqs. 03]) and (gj) can be written 
as 



SF" 



X3 ~ 7^1 
fc 2 „ 

4 



SL 



Xx-X2 - x 3 )S v 



3xi 



3zi - 3xi(x 3 - 1) >SF 



+3(-2xi + x 3 - 1)5F' 
+ (1 - 2 Xl + x 3 )5F" 



35F' 1 



(63) 



12 



rz-2 



^■^3 / / 2 

2 — 2x 3 H Xl(X3 + 1) — x 1 + x 1 



= (1 — xi - x 2 - x 3 )S m - xiS' m , 
where a new variable, £5 = aH , satisfies 

x' 5 = (x 3 - l)x 5 . 
Note that the perturbation SR is given by 

6R=6H 2 ^6F. 
m 



SF 



(64) 



(65) 



(66) 



The evolution of the Hubble parameter is known by solv- 
ing the equation 



H 



X3 



(67) 



Solving equations (|63|) and (|64[1 together with the back- 
ground equations (|2"5|) - (f50)) . (p5|) and (f6"T)) numerically, 
we find the evolution of 8 rn and <5i?. 
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